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TECHNICAL  SUMMARY 


The  objective  of  this  project  is  to  determine  the  yield  bias  of  underground  nu¬ 
clear  tests  induced  by  the  presence  of  a  high  velocity  descending  slab  beneath 
the  test  site.  Specifically,  the  effect  of  the  Aleutian  slab  is  being  investigated 
on  the  US  underground  tests  Longshot,  Milrow,  and  Cannikan.  P  wave  seis¬ 
mograms  will  be  synthesized  using  dynamic  ray  tracing  and  superposition 
of  Gaussian  beams  in  three-dimensional  models  of  the  Aleutian  slab  deter¬ 
mined  from  P  travel  time  delays.  Focusing  and  defocusing  and  multipathing 
at  teleseismic  distances  will  be  evaluated  by  comparison  of  observed  with 
synthetic  seismograms  of  the  Aleutian  tests. 

Data  collection  for  of  Amchitka  P  waveforms  and  amplitudes  was  ini¬ 
tiated.  A  no  cost  extension  of  this  project  was  requested  and  granted  by 
AFGL  to  allow  for  collection  of  the  necessary  data  to  complete  the  project. 
Pending  completion  of  the  data  collection,  work  was  begun  on  a  project  to 
compute  complete  regional  seismograms  in  crustal  and  upper  mantle  models 
having  gradients  in  layers.  This  project  will  be  included  in  a  new  proposal 
and  directly  contributes  to  the  AFGL  treaty  verification  program,  with  its 
current  emphasis  on  CTBT  monitoring  at  local  and  regional  distances. 

In  collaboration  with  Danny  Harvey,  the  locked  mode  method  of  syn¬ 
thesizing  complete  regional  seismograms  (Harvey,  1981)  was  modified  to  in¬ 
clude  the  Langer  uniform  asymptotic  approximation  to  vertical  wavefunc- 
tions  within  layers  having  linear  vertical  velocity  gradients.  Good  agree¬ 
ment  is  obtained  in  gradient  models  between  synthetics  computed  using  the 
Langer-locked  mode  method,  the  colocation  method,  and  the  conventional 
locked  mode  method  in  models  parameterized  by  thin  homogeneous  layers. 
Errors  in  calculated  displacement  introduced  by  the  use  of  the  Langer  ap¬ 
proximation  remain  less  than  several  percent  for  wavelengths  A  <  0.2 V/W. 
Whenever  it  is  necessary  to  represent  gradients  accurately,  the  Langer-locked 
mode  method  is  computationally  more  efficient  than  the  locked  mode  method 
using  thin  homogeneous  layers.  By  reducing  the  number  of  parameters 
needed  to  describe  an  Earth  model,  the  Langer-locked  mode  method  will 
also  simplify  the  inverse  problem  of  determining  structure  using  observed 
and  synthetic  regional  seismograms.  Test  calculations  of  regional  seismo¬ 
grams  confirm  that  the  Pn  and  Sn  phases  are  strongly  affected  by  the  mag¬ 
nitude  of  the  velocity  gradients  beneath  the  Moho,  but  that  Lg  is  only  weakly 
*fiected  by  the  details  of  crustal  layering. 
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SYNTHESIS  OF  COMPLETE  SEISMOGRAMS  BY  THE  LOCKED  MODE  METHOD 
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ABSTRACT 

Any  realistic  crustal  and  upper  mantle  model  possesses  layers  with  vertical  gradients.  Elastic  mod¬ 
uli  and  density  in  each  layer  are  affected  by  pressure,  temperature,  pore  fluids,  and  crack  density. 
All  of  these  quantities  change  continuously  with  depth,  many  having  a  well  known  functional  depen¬ 
dence  on  depth.  Virtually  all  of  the  regional  phases  can  be  strongly  affected  by  velocity  gradients. 
The  best  known  effects  of  velocity  gradients  are  on  the  Pn  and  Sn,  in  which  small  changes  in  the 
velocity  gradient  beneath  the  Moho  can  make  large  changes  in  the  decay  of  Pn  and  Sn  with  dis¬ 
tance.  Methods  of  synthesizing  complete  regional  seismograms  often  inadvertently  ignore  the  effect 
of  crustal  gradients  by  parameterizing  the  Earth  model  with  planar  homogeneous  layers.  To  remedy 
this  problem  we  have  modified  the  locked  mode  method  of  synthesizing  complete  regional  seismo¬ 
grams  to  include  the  Langer  uniform  asymptotic  approximation  to  vertical  wavefunctions  within 


layers  having  linear  vertical  velocity  gradients.  Good  agreement  is  obtained  in  gradient  models 
between  synthetics  computed  using  the  Langer-locked  mode  method,  the  colocation  method,  and 
the  conventional  locked  mode  method  in  models  parameterized  by  thin  homogeneous  layers.  Errors 
in  calculated  displacement  introduced  by  the  use  of  the  Langer  approximation  remain  less  than 
several  percent  for  wavelengths  A  <  Q.2V/W.  Whenever  it  is  necessary  to  represent  gradients 
accurately,  the  Langer-locked  mode  method  is  computationally  more  efficient  than  the  locked  mode 
method  using  thin  homogeneous  layers.  By  reducing  the  number  of  parameters  needed  to  describe 
an  Earth  model,  the  Langer-locked  mode  method  will  also  simplify  the  inverse  problem  of  deter¬ 
mining  structure  using  observed  and  synthetic  regional  seismograms.  Test  calculations  of  regional 
seismograms  confirm  that  the  Pn  and  Sn  phases  are  strongly  affected  by  the  magnitude  of  the 
velocity  gradients  in  beneath  the  Moho,  but  that  Lg  is  only  weakly  affected  by  the  details  of  crustal 
layering. 


INTRODUCTION 

Complete  seismograms  at  local  and  regional  distances  are  now  routinely  computed  in  plane  layered 
models  for  a  variety  of  source  receiver  geometries,  source  depths,  and  source  types  by  integrating 
or  summing  over  wavenumbers  (Bouchon  and  Aki,  1977;  Kind  1978;  Wang  and  Hermann,  1980) 
or  summing  locked  or  leaky  modes  (Harvey,  1981;  Kerry,  1981;  Haddon,  1986;  Nolet  et  al.,  1989). 
The  computational  expense  of  these  calculations  remains  relatively  cheap  as  long  as  the  crust  and 
upper  mantle  model  can  be  described  by  a  small  number  of  planar,  homogeneous  layers. 

Seismograms  synthesized  in  models  composed  of  small  number  of  plane  homogeneous  layers 
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ignore  the  continuous  depth  dependence  of  elastic  moduli.  Usually  seismograms  are  synthesized  in 
simple  models  composed  of  a  two  or  three  homogeneous  layers  of  crust  overlaying  a  homogeneous 
lid,  low  velocity  zone,  and  upper  mantle  beneath  the  lid.  Since  Earth  curvature  is  ignored  in  these 
calculations,  the  model  is  effectively  one  in  which  each  layer  has  a  small  negative  gradient  with 
depth. 

The  simplest  generalization  of  a  homogeneously  layered  model  is  to  allow  for  the  effect  of  velocity 
gradients.  Any  realistic  crustal  and  upper  mantle  model  possesses  layers  with  vertical  gradients. 
Elastic  moduli  and  density  in  each  layer  are  affected  by  pressure,  temperature,  pore  fluids,  crack 
density  and  aspect  ratio.  All  of  these  quantities  change  continuously  with  depth,  many  having  a 
well  known  functional  dependence  on  depth. 

Virtually  all  of  the  regional  phases  can  be  strongly  affected  by  velocity  gradients.  One  example 
is  dispersion  of  the  fundamental  mode  Rayleigh  wave,  or  Rg  phase  at  local  and  regional  distances. 
Crustal  models  having  a  homogeneous  layer  at  the  surface  produce  an  unrealistically  impulsive, 
undispersed  Rg  arrival.  To  match  observed  data,  the  fundamental  mode  arrival  must  be  artificially 
removed  or  attenuated.  Perhaps  the  best  known  effects  of  gradients  on  regional  phases  are  those 
on  the  Pn  and  Sn  phases.  In  a  plane  layered  model,  the  Pn  and  Sn  phases  are  classical  head  waves 
traveling  just  beneath  the  Moho.  Hill  (1971)  and  Cerveny  and  Ravindra  (1971)  have  shown  how 
gradients  transform  classical  head  waves  into  interference  head  waves  or  ’’whispering  gallery  waves” 
(e.g.,  Cormier  and  Richards,  1976;  Menke  and  Richards,  1980).  The  distance  decay  of  both  classical 
and  interference  headwaves  is  frequency  dependent. 

In  this  paper,  we  describe  the  results  of  incorporating  velocity  gradients  in  crustal  and  upper 


mantle  models  using  the  locked  mode  method.  Gradients  are  introduced  into  the  computations 
by  allowing  each  layer  to  be  vertically  inhomogeneous  and  applying  the  Langer  approximation 
(Appendix  I)  to  calculate  an  asymptotic  approximation  in  frequency  to  the  vertical  wavefunctions 
in  each  inhomogeneous  layer. 

A  brief  review  of  the  locked  mode  method  is  first  given.  Mathematical  details  of  the  Langer 
approximation  and  its  incorporation  in  the  locked  mode  method  are  described  in  Appendices  I  and 
II.  The  remainder  of  the  paper  describes  the  results  of  tests  conducted  to  determine  the  accuracy 
of  the  Langer  approximation  and  how  it  breaks  down  as  the  gradient  in  the  layer  increases.  A 
discussion  and  example  show  how  depth  and  frequency  dependent  attenuation  can  be  included  in 
the  Langer-locked  mode  method.  The  paper  concludes  with  discussion  of  synthetic  seismographs 
showing  how  gradients  near  the  free  surface  and  Moho  can  radically  affect  the  propagation  of  some 
of  the  principal  regional  phases. 

Review  of  the  Locked  Mode  Method 
Following  Harvey  (1981;  1985),  the  complex  displacement  spectra  are  evaluated  from 
Ru(u>,xT,0r,zr)  ~  -RvI-RpI-iYlYl  flA(rc.w)  RET{n,u,m)  RE  (n,u,z,)  R\p(n,u;,m,  xr,9T,  zr) 

n  m 

(i) 

Lu(u,xr,0r,Zr)  =  ~LpI  -  'J2J2  £A(n-w)  T(n>  m)  LE(n,w,r,)  Lilin,u,m,xT,6r,Zr) 

n  m 

where  the  subscripts  R  and  L  denote  Rayleigh  and  Love  modes  respectively. 
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rA  and  IjA  area  scalar  amplitude  factors  defined  by 


flA(n,u>)  =  - 


**y23(o) 

dYu(Q)/dk 


(2) 


lA{ti,lj)  =  - 


tkD2(  0) 


dDi(0)/dk 

rT,t  and  £,Et  are  row  vectors  defined  from  the  source  jump  vectors,  r tp  and  i tj)  are  defined  from 
products  of  eigenfunctions  for  displacement  (rE\,  rE2)  and  lE)\  and  vector  cylindrical  harmonics 
P,B,  C: 


RTp(n,u,m,xr,0r,zr)  =  REi(n,bj,zr)P(n,u>,m,xr,0r)  +  RE2(n,u,  zr)B(n,u,m,xr,0r) 


(3) 


Li){n,u,m,xT,0r,zT )  =  LEi(n,u>,zT)C(n,u;,m,xr,0T ) 


r^I  and  LgI  are  branch  cut  integrals,  which  account  for  energy  that  cannot  be  represented  by 
normal  modes,  and  are  associated  with  near  vertically  propagating  P  and  S  waves  that  leak  into 
the  halfspace.  The  locked  mode  method  does  not  evaluate  the  branch  cut  integrals.  It  chooses 
the  halfspace  to  be  sufficiently  deep  and  fast  such  that  all  of  the  energy  important  to  a  particular 
time  window  at  a  particular  distance  can  be  accurately  represented  by  the  locked  mode  summation 
alone. 

Seismograms  are  synthesized  by  evaluating  the  complex  spectra  at  discrete  frequencies  and 
inverting  to  the  time  domain  by  fast  Fourier  transform.  A  small  complex  frequency  can  be  added 
to  attenuate  all  arrivals  that  arrive  outside  of  the  finite  time  window  given  by  the  folding  frequency 
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of  the  discrete  Fourier  transform  (Rosenbaum,  1974;  Muller  and  Schott,  1981).  Harvey  (1981, 
1985)  gives  detailed  derivations  of  the  locked  mode  method  and  describes  its  implementation  in 
media  described  by  homogeneous  layers.  The  principal  modifications  of  the  method  for  use  with 
the  Langer  approximation  are  concerned  with  the  calculation  of  the  eigenfunction  vector  E  and 
the  scalar  amplitude  factors  /?A  and  iA.  The  partial  derivatives  with  respect  to  k  appearing  in 
the  amplitude  factors  are  calculated  by  difference  derivatives.  Appendices  I  and  II  describe  the 
calculation  with  the  Langer  approximation  of  the  Y  matrix  elements  and  the  vectors  D  and  E. 

The  Langer  approximation  can  also  be  implemented  in  methods  of  synthesizing  complete  seis¬ 
mograms  that  numerically  integrate  over  horizontal  wavenumber  and  slowness  (Cormier,  1980). 
The  primary  advantage  of  the  locked  mode  method  is  that  most  of  the  computational  effort  in¬ 
volved  in  the  calculation  of  the  amplitude  factors  and  eigenfunctions  can  be  catalogued  for  use  with 
different  source-receiver  geometries  and  different  moment  tensor  representations  of  point  sources. 
Although  response  functions  can  be  similarly  cataloged  in  approaches  that  integrate  or  sum  over 
wavenumber  or  slowness,  this  is  rarely  done  in  practice.  A  secondary  advantage  of  the  locked  mode 
method  is  that  a  large  body  of  literature  exists  in  modal  notation  on  inversion  for  structure  and 
source  parameters.  The  analysis  of  problems  using  normal  modes  of  the  whole  Earth  at  low  fre¬ 
quency  and  long  range  can  usually  be  directly  adapted  to  higher  frequency  and  shorter  range  using 
locked  modes  (e.g.,  Gomberg  and  Masters,  1988). 
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The  Accuracy  of  the  Langer  Approximation 


The  Lange*  approximation  assumes  decoupling  between  P  and  S  waves  and  up-  and  down-going 
waves  in  each  gradient  layer,  and  the  criteria  for  its  accuracy  are  thus  similar  to  those  used 
in  ray-asymptotic  solutions  to  the  elastodynamic  equation  of  motion  in  inhomogeneous  media 
(Richards,  1976).  Qualitatively,  the  Langer  approximation  is  known  to  become  less  accurate  as 
non-dimensional  ratios  A/(n/Vu)  increase,  where  v  is  a  P  or  S  velocity  or  density  (Richards,  1976; 
Chapman,  1974).  Another  way  in  which  this  is  commonly  phrased  is  that  the  wavelength  must  be 
much  smaller  than  the  scalelength  of  the  medium,  /,  where  l  is  the  maximum  of  (a/Va,  /3/V/3,p/Vp) 
(Beydoun  and  Ben-Menahem,  1985).  A  goal  in  this  study  was  to  quantify  the  breakdown  in  the 
Langer  approximation  as  the  scalelength  of  gradient  layers  decrease,  determining  exactly  how  large 
the  ratio  X/l  can  be  before  errors  in  calculated  displacement  exceed  some  specified  bound. 

The  first  step  in  such  a  study  is  to  choose  accurate  reference  synthetic  seismograms  in  models 
having  strong  gradients.  Spudich  and  Ascher  (1983)  published  synthetic  seismograms  calculated 
by  the  numerical  colocation  method  for  a  simple  model  consisting  of  a  gradient  over  half  space. 
The  gradient  layer  in  this  model  was  parameterized  by  a  sequence  of  40  thin  layers  (Figure  1),  the 
width  of  each  thin  layer  approximately  equal  to  one-tenth  the  wavelength  of  shear  waves  at  1  Hz. 
Excellent  agreement  was  found  between  the  locked  mode  synthetics  and  the  colocation  synthetics. 
This  result  confirmed  that  locked  mode  synthetics  computed  in  models  in  which  gradient  layers 
are  represented  by  thin  layers  can  be  used  as  accurate  reference  synthetics  to  test  the  Langer 
approximation. 

To  test  the  accuracy  of  the  Langer  approximation,  seismograms  were  synthesized  using  the 
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locked  mode  method  using  the  Langer  approximation  in  a  series  of  models  with  increasing  gradients 
in  P  and  S  velocity  and  density  in  a  layer  over  a  halfspace  (Figure  2).  Figure  3  compares  the 
dispersion  curves  of  the  locked  Love  and  Rayleigh  modes  calculated  with  Langer  approximation 
in  a  thick  continuous  gradient  layer  with  those  calculated  by  parameterizing  the  gradient  layer 
with  thin  homogeneous  layers.  Even  for  the  mo6t  severe  of  the  gradients  shown  in  Figure  2,  the 
dispersion  curves  calculated  using  the  Langer  approximation  remain  quite  accurate  throughout 
nearly  the  entire  range  of  phase  velocity  and  frequency.  The  primary  region  of  error  occurs  for  the 
low  frequencies  of  the  fundamental  mode.  This  is  not  unexpected  since  most  of  the  energy  of  the 
fundamental  mode  in  this  frequency  band  is  confined  to  the  strong  gradient  layer  near  the  surface. 
As  expected,  the  errors  in  the  dispersion  curves  calculated  by  the  Langer  approximation  are  largest 
at  low  frequency,  where  the  wavelength  approaches  the  scale  length  of  the  gradient  layer. 

Figure  4  compares  reference  synthetics  and  Langer  approximated  synthetics  for  the  sequence  of 
gradient  models  shown  in  Figure  2.  Although  the  kinematic  errors  in  the  mode  dispersion  calcula¬ 
tions  are  small  throughout  most  of  the  frequency  band,  the  dynamic  errors  in  mode  amplitudes  are 
sufficient  to  produce  poor  matches  in  the  group  velocity  band  corresponding  to  the  fundamental 
mode  and  the  first  few  higher  modes.  These  effects  can  be  seen  in  Figure  4,  in  which  the  earlv 
portion  of  the  seismograms  computed  by  the  two  methods  are  more  closely  in  phase  but  become 
progressively  out  of  phase  in  the  time  window  corresponding  to  the  arrival  of  the  fundamental 
mode  and  first  few  higher  modes.  The  agreement  between  the  two  methods  is  much  better  for  the 
transverse  component  than  the  radial  or  vertical  components  of  motion. 

The  match  between  reference  and  Langer  approximated  synthetics  becomes  nearly  perfect  for 


weakest  surface  gradients  (model  3  in  Figure  2).  The  seismograms  computed  by  the  two  methods 
overlay  one  another  to  within  the  thickness  of  plotted  lines.  The  difference  seismograms  in  Figure 
5  are  largest  near  the  peak  oscillations  where  small  differences  in  arrival  time  of  pulses  having  high 
slopes  produce  large  differences.  Since  the  dominant  frequency  the  synthetic  seismograms  is  about 
0.5  Hz.,  one  can  conclude  that  errors  in  the  use  of  the  Langer  approximation  become  less  than 
several  percent  when  the  ratio  A//  is  less  than  or  equal  to  0.2.  If  one  were  not  interested  in  the 
accuracy  of  the  fundamental  mode  at  frequencies  less  than  1  Hz,  the  Langer  approximation  could 
synthesize  the  higher  modes  in  this  example  with  high  accuracy  across  the  entire  frequency  band. 
The  fundamental  mode  could  be  synthesized  with  high  accuracy  at  frequencies  greater  than  1  Hz. 

Intrinsic  Attenuation 

To  be  practically  useful,  any  method  of  synthesizing  complete  seismograms  at  local  and  regional 
distances  must  be  capable  of  including  intrinsic  attenuation.  The  incorporation  of  the  attenuation 
in  the  Langer  approximation  simply  consists  of  the  analytic  continuation  of  all  formulae  to  complex 
velocities  (Cormier  and  Richards,  1976,  1989).  Care  must  be  exercised  in  the  definition  of  branch 
cuts  of  square  roots  and  fractional  powers  appearing  in  both  the  analytic  expressions  and  function 
subroutines  used  in  evaluating  the  Langer  approximation  (see  Appendix  I),  but  this  is  not  an 
insurmountable  problem.  The  Langer  subroutine  modified  for  use  with  locked  mode  calculations 
has  been  tested  in  problems  involving  integration  in  the  complex  ray  parameter  plane  combined 
with  complex,  frequency  dependent  velocity.  It  returns  generalized  vertical  wavefunctions  and 
slownesses  that  are  continuous  in  the  complex  ray  parameter  plane  except  for  poles  and  branch 
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cuts,  which  emanate  from  complex  ray  parameters  corresponding  to  grazing  incidence  on  boundaries 
in  an  anelastic  model.  Test  calculations  have  demonstrated  that  the  position  of  these  singularities 
do  not  impede  a  successful  search  for  the  complex  zeros  of  the  dispersion  functions  of  locked  modes 


where 


a  given  region.  Complex  velocities  are  calculated  at  each  layer  boundary  by  equation  5  above 
and  linear  gradients  of  complex  velocity  are  are  assumed  in  each  layer.  The  delay  time  function  r 
needed  by  the  Langer  approximation  is  calculated  as  described  in  Appendix  II,  but  it  now  must 
be  recalculated  at  each  frequency.  It  is  possible  to  specify  different  peak  Qp  values  as  well  as 
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different  upper  and  lower  limits,  ui\  and  0J2,  of  the  relaxation  band  at  the  top  and  bottom  of  each 
inhomogeneous  layer. 

A  test  anelastic  model  is  shown  in  Figure  6.  The  attenuation  model  is  an  absorption  band  model 
in  pure  shear  attenuation  having  gradients  in  peak  attenuation  Qp1,  and  low  and  high  frequency 
corners,  wj.wj,  of  the  relaxation  band.  A  minimum  value  of  Qp  =  20  is  assumed  at  the  surface. 
The  velocities  and  Q  values  are  similar  to  values  measured  from  regional  seismograms  in  New 
England  (Kafka  and  Reiter,  1987).  Locked  mode  seismograms  were  synthesized  in  these  model  using 
two  different  approaches.  In  the  first  approach,  only  the  real  part  of  the  complex  velocities  was  used 
in  calculation  of  mode  amplitudes  and  eigenfunctions,  a  complex  phase  velocity  was  substituted  in 
the  cylindrical  harmonics  describing  the  horizontal  propagation  of  each  mode.  This  complex  phase 
velocity  is  taken  from  the  complex  pole  k  estimated  by  first  order  perturbation  theory.  This  is  the 
standard  approach  for  handling  attenuation  in  surface  wave  and  locked  mode  calculations  (Harvey, 
1985;  Panza  and  Sudhadolc,  1987),  and  is  assumed  to  be  accurate  if  the  Q  factor  is  sufficiently  high. 
Day  et  al.  (1989)  have  shown  this  approach  to  be  inaccurate  for  some  regional  seismic  phases  even 
at  Q  values  on  the  order  of  several  hundred.  For  this  reason,  seismograms  were  also  synthesized  by 
an  exact  approach,  in  which  a  search  was  made  for  the  complex  roots  of  the  dispersion  function  and 
all  formulae,  including  ampiitude'factors  and  eigenfunctions,  were  evaluated  at  these  complex  roots. 
The  complex  pole  searching  algorithm  was  based  on  one  suggested  by  Schwab  and  Knopoff  (1971), 
with  modifications  near  osculating  points  of  the  dispersion  curves.  Near  these  points,  the  complex 
roots  are  found  by  the  same  algorithm  for  a  series  of  increasing  Q~x  values,  approaching  the  true 
Q~l  model.  Checks  are  made  for  duplication  or  omission  of  poles  at  the  end  of  this  procedure  for 


11 


each  frequency. 

Figure  7  compares  the  results  of  these  two  methods  for  incorporating  attenuation  of  the  funda¬ 
mental  mode  Rayleigh  wave.  The  seismograms  computed  by  the  different  methods  nearly  overlay 
one  another  at  all  distances.  The  exact  method  reduces  some  high  frequency  numerical  noise, 
which  is  barely  visible  at  the  scale  of  Figure  7.  The  differences  in  the  complex  phase  velocities 
computed  by  the  two  methods  are  on  the  order  of  0.001  km/sec  in  the  real  part  of  the  complex 
phase  velocity  and  vary  from  1  x  10~10  to  1  *  10“4  km/sec  in  the  complex  part  of  phase  velocity  as 
frequency  increases  up  to  2  Hz.  The  differences  between  the  depth  behavior  of  the  real  part  of  the 
complex  eigenfunctions  are  insignificant  between  the  two  methods.  From  these  results  it  can  be 
concluded  that  the  perturbation  approach  to  attenuation  remains  very  accurate  in  the  synthesis  of 
the  fundamental  mode  for  Q  values  as  low  as  20.  For  the  synthesis  of  higher  modes,  particularly 
those  contributing  to  refracted  P  and  S  and  interference  head  waves,  more  detailed  tests  have  shown 
that  the  perturbation  approach  introduces  significant  error  as  Q  values  decrease  below  several  100. 

It  is  reasonable  to  assume  that  gradients  in  the  real  part  of  elastic  moduli  are  also  associated 
with  gradients  in  the  imaginary  part  of  elastic  moduli.  We  have  demonstrated  in  this  section  that 
the  Langer  approximation  can  be  applied  to  locked  mode  calculations  in  models  having  gradients 
in  complex  elastic  moduli.  Often  a  very  low  Q  layer  is  required  in  a  surface  layer  in  order  to 
produce  realistic  simulations  of  seismograms  observed  at  local  and  regional  distances  (e.g.,  Panza 
and  Sudhadolc,  1987).  If  the  apparent  attenuation  of  such  a  layer  is  truly  due  to  viscoelasticty, 
its  effects  can  be  accurately  calculated  by  complex  locked  modes.  It  is  worth  noting,  however, 
that  such  apparent  low  Q' s  are  likely  due  to  a  combination  of  scattering  by  topography  of  layer 
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boundaries  and  volumetric  heterogeneities  and  frictional  sliding  of  grains  and  open  cracks.  Neither 
of  these  effects  can  be  simulated  by  a  combination  of  vertically  varying  layers  and  linear  viscoelastic 
relaxations. 


Effects  of  Gradients  in  Realistic  Models 

To  test  the  effects  of  crustal  and  upper  mantle  gradients  on  regional  seismic  phases,  locked  mode 
synthetics  were  computed  in  two  simple  models  MH  and  MG  (Figure  8).  Model  MH  consists  of 
a  two-layered  crust  overlaying  a  homogeneous  mantle.  MH  has  also  been  used  for  testing  and 
benchmark  timing  of  many  different  techniques  of  computing  complete  seismograms  at  regional 
distances  (Richards  and  Mithal,  personal  communications)  Model  MG  consists  of  a  single  crustal 
layer  having  a  positive  gradient  with  depth,  overlaying  a  mantle  having  a  positive  gradient  with 
depth.  The  mantle  gradient  is  consistent  with  the  increase  in  seismic  velocities  typical  of  reference 
earth  models  between  the  Moho  and  400  km  depth.  The  depth  averaged  crustal  velocities  of  MH 
and  MG  are  identical.  Both  models  have  an  attenuation  structure,  with  Q’s  in  a  high  enough  range 
that  simple  perturbation  theory  can  be  accurately  used  to  calculate  the  effects  of  attenuation  in 
the  locked  mode  method.  Seismograms  were  synthesized  in  a  frequency  band  up  to  2  Hz.  for  the 
source  and  receiver  geometries  used  by  W-Y.  Kim  (1987),  who  synthesized  seismograms  in  model 
MH  using  wavenumber  integration. 

The  synthetic  seismograms  for  the  first  10  Rayleigh  modes  (Figure  9)  are  very  similar  for  both 
model  MH  and  MG.  The  group  velocity  window  of  the  energy  centroid  corresponds  to  that  expected 
for  the  Lg  phase.  The  strong  similarity  of  the  synthetic  seismograms  suggests  that  Lg  is  not  very 


13 


sensitive  to  the  details  of  the  crustal  model,  its  coda  primarily  being  controlled  by  the  total  thickness 
of  the  crust  and  its  average  shear  velocity.  It  is  probably  possible  to  simulate  realistic  Lg  phases 
using  a  very  few  number  of  crustal  layers.  Introduction  of  crustal  layers  in  a  modeling  experiment 
may  not  be  necessary  unless  there  is  compelling  evidence  for  crustal  discontinuities  observed  in  the 
earlier  time  window  in  the  form  of  refracted  body  waves  and  interference  head  waves. 

In  a  comparison  of  complete  seismograms  (Figure  10),  the  seismograms  are  very  similar  at 
closer  ranges  but  at  300  km  some  differences  begin  to  be  notable.  Pn  and  Sn  are  very  weak  in 
the  MH  simulation,  but  are  very  strong  in  the  MG  synthetic.  Pn,  Sn,  and  crustal  reverberations 
converted  to  Pn  and  Sn  are  so  strong  in  the  MG  synthetic  that  they  dominate  Lg  in  amplitude, 
and  the  seismogram  seems  to  be  a  series  of  spikes  when  the  display  is  scaled  on  the  peak  amplitude. 
The  comparison  confirms  what  is  known  about  Pn  and  Sn  as  interference  head  waves  in  models 
having  positive  gradients  below  the  Moho.  A  positive  gradient  acts  to  enhance  the  amplitude 
of  the  interference  head  wave  far  above  what  would  be  predicted  for  a  classical  head  wave  in  a 
homogeneous  layer  (Hill,  1971;  Cerveny  and  Ravindra,  1971;  Menke  and  Richards,  1980). 

CONCLUSIONS 

In  this  paper,  we  have  demonstrated  that  even  small  gradients  of  VK  =  0.03  sec.-1  can  substantially 
affect  the  distance  decay  of  interference  head  waves  such  as  Pn  and  Sn.  Lg,  on  the  other  hand,  is 
only  verly  weakly  sensitive  to  details  of  crustal  layering  or  gradients.  The  peak  amplitude  and  coda 
length  of  Lg  primarily  depends  on  total  crustal  thickness  and  average  shear  velocity  of  the  crust. 

In  either  a  locked  mode  or  wavenumber  integration  approach  to  synthesizing  complete  seismo- 
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grams,  the  Langer  approximation  can  accurately  approximate  vertical  wavefunctions  in  inhomoge¬ 
neous  layers  having  a  single  ray  turning  point  for  wavelengths  that  are  small  with  respect  to  the 
scale  length  of  the  layer.  This  result  can  be  quantified  by  stating  that  errors  in  the  amplitude  and 
phase  of  synthetic  seismograms  are  less  than  several  per  cent  for  wavelengths  A  <  0.2V/VV,  where 
Via  velocity  or  density  function.  At  5  Hz.  this  inequality  is  statisfied  by  gradients  beneath  the 
Moho  as  high  as  0.8  sec.-1. 

Propagation  of  the  wavefield  using  the  Langer  approximation  in  a  vertically  inhomogeneous 
layer  will  often  represent  a  computational  savings  over  propagation  though  the  gradient  layer  pa¬ 
rameterized  by  a  sequence  of  thin  layers.  An  example  of  a  gradient  layer  parameterized  by  40  thin 
homogeneous  layers  executed  about  a  factor  of  two  slower  in  both  the  pole  searching  and  eigenfunc¬ 
tion  evaluation  compared  to  the  same  calculations  using  the  Langer  approximation  in  the  gradient 
layer  parameterized  by  analytic  velocity  functions.  A  calculation  in  a  thick  homogeneous  layer, 
however,  would  still  always  be  more  efficient  than  a  calculation  using  the  Langer  approximation  in 
an  inhomogeneous  layer  of  the  same  thickness.  A  model  parameterization  that  may  be  the  best 
compromise  between  computational  efficiency  and  realism  in  the  behavior  of  regional  phases  would 
be  one  having  a  crust  composed  of  homogeneous  layers  overlaying  a  mantle  composed  of  gradient 
layers.  Seismograms  synthesized  in  such  a  model  could  accurately  predict  the  Lg  phase  as  well  as 
the  Pn  and  Sn  phases.  (This  study  did  not  investigate  the  importance  of  crustal  gradients  for  the 
Pg  phase.) 

The  Langer  locked  mode  approach  to  synthesizing  complete  seismograms  may  also  offer  some 
advantages  in  waveform  inversion  for  earth  structure.  By  reducing  the  number  of  parameters  needed 
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to  describe  a  model,  the  inverse  problem  for  structure  would  be  simplified  and  fewer  experiments 
would  be  needed  to  determine  the  maximum  number  of  resolvable  layers.  A  layer  need  only  be 
introduced  whenever  the  data  firmly  suggest  the  existence  of  first  order  discontinuities. 
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APPENDIX  I  -  THE  LANGER  APPROXIMATION 


J 

Vertical  Wavefunctions  j 


The  notation  for  the  Langer  approximation  (Langer,  1932;  1949)  differs  among  different  authors 
who  have  applied  it  to  seismic  wave  propagation.  (Richards,  1976;  Woodhouse,  1978;  Chapman, 
1974;  Doornbos,  1981),  involving  either  Hankel  functions  of  order  1/3  or  Airy  functions  of  different 


types  or  arguments  to  give  exponentially  decaying  and  growing  type  solutions  below  a  turning 
point.  The  notation  adopted  here  is  basically  that  given  in  Aki  and  Richards  (1980). 

The  Langer  approximation  is  a  uniformly  asymptotic  approximation  to  the  vertically  separated 
part  of  the  solution  to  the  elastodynamic  wave  equation  in  a  region  in  which  elastic  moduli  and 
density  vary  continuously  with  depth.  The  zeroth  order  term  in  frequency  in  the  asymptotic  solution 


is  given  as 


7T^(r)  =  V2  7 r  e  3 

fey41 

uyy 

j  Ai(—  e  3 
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where 


Ai  is  an  Airy  function  and 


C*  =  (3/2ura)V* 


C/3  ==  (3/2  urp)2'3 


XQ  ~ 


Xp  =  y/lfl3'-  ?fr* 


a  and  /?  are  the  P  and  S  velocity  respectively  at  radius  r,  p  is  the  ray  parameter  in  a  spherical 
Earth,  and  rp  is  the  turning  point  radius,  i.e.,  that  radius  at  which  Xa  or  Xp  vanishes.  In  each 
inhomogeneous  layer,  the  velocity  functions  o(r)  and  (3{t)  are  assumed  to  be  analytic  and  to 
possess  only  one  turning  point  rp  in  the  domain  of  complex  p  used  in  synthesizing  a  seismogram. 
The  Langer  approximated  wavefunctions  can  also  be  written  in  terms  of  Hankel  functions  of  order 
1/3  (Richards,  1976;  Doornbos,  1981). 

The  ir  wavefunctions  are  those  for  P  waves;  the  a  wavefunctions  are  those  for  S  waves.  Several 
possible  pairs  of  independent  solutions  may  be  chosen  to  define  fundamental  matrices,  which  can 
be  used  to  solve  problems  in  wave  propagation  in  media  consisting  of  a  sequence  of  vertically 
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inhomogeneous  layers.  The  pairs  (x^,  t(2))  and  <r<2))  correspond  to  up-  (1)  and  down¬ 

going  (2)  waves.  The  pairs  (x*1),  x<3))  and  (ffM,  cr(3))  correspond  to  up-going  (1)  and  standing  or 
evanescent  waves  (3).  When  the  turning  point  radius  rp  is  greater  than  r,  the  wavefunctions  x<3> 
and  o^3)  are  always  exponentially  decaying  functions  with  decreasing  radius  r 


Vertical  Slownesses 


Implementation  of  the  Langer  approximation  in  problems  in  which  elastic  boundary  conditions 
must  to  be  satisfied  at  model  discontinuities  is  simplified  by  the  introduction  of  generalized  cosines 
(Richards,  1976;  Aki  and  Richards,  1980)  or  generalized  vertical  slowness  functions,  which  are 
defined  as  follows 


dxM 

dr 


/(iW1*) 


* 

* 

z 


^r/(,W7r(1)) 

— —  /(tunW) 


V  = 

V  = 

V  = 


doW 

dr 

d<jW 

dr 

daW 

dr 


(AI.2) 


The  normalization  of  the  vertical  wavefunctions  differs  slightly  from  that  given  in  Aki  and 
Richards  (1980)  and  has  been  chosen  such  that  the  following  relations  are  satisfied 


6r(1M2> -KW1)  =  1 
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^(Dt(3)+^(3)t(1)  =  i 


(AI.3) 

+  r,^)^)  =  1 

+  fja^x^  =  1 

These  relations  can  be  demonstrated  by  substituting  the  Langer  approximation  to  the  vertical 
slownesses  and  the  Wronskian  relations  between  the  Airy  functions  having  different  arguments. 
Equations  AI.3  are  satisfied  exactly  when  only  the  zero  order  terms  in  frequency  are  kept  in  the 
definitions  of  the  vertical  slownesses. 

Branch  Cuts  and  Complex  p 

The  functions  that  define  the  generalized  vertical  wavefunctions  and  slownesses  as  well  as  the  special 
function  subroutines  from  which  Airy  functions  or  Hankel  functions  of  order  1/3  are  commonly 
constructed  contain  branch  cuts  emanating  from  points  in  the  complex  p  plane  corresponding 
to  ray  parameters  grazing  the  model  discontinuities.  Extreme  care  must  be  exercised  both  in  the 
definition  and  the  choosing  of  branch  cuts  appearing  in  all  functions  of  variables  raised  to  fractional 
powers.  A  subroutine  for  the  Langer  approximated  wavefunctions  and  vertical  slownesses  has  been 
used,  in  which  the  branch  cuts  of  the  vertical  wavefunctions  are  defined  as  in  Cormier  (1976)  and 
Cormier  and  Richards  (1989).  This  wavefunction  subroutine  has  been  tested  in  a  wide  variety 
of  problems  involving  both  complex  velocities  and  complex  p.  For  examples,  see  Cormier  and 
Richards,  (1989). 
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Fundamental  Matrices 


Boundary  conditions  in  an  medium  consisting  of  n  inhomogeneous  layers  can  be  handled  in  the 
same  manner  as  a  medium  consisting  of  homogeneous  layers,  but  with  the  Langer  approximation 
to  the  vertical  wavefunctions  and  vertical  slownesses  substituting  for  exponential  functions  and 
cosines. 


P-SV 


As  a  function  of  radius  r,  the  fundamental  matrix  for  P-SV  propagation  and  Rayleigh  modes  is 
taken  to  be  than  given  in  Cormier  (1980): 
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— iBfjo ^  - Au W  p/rer(l)  irjcr^  j 

The  fundamental  matrix  may  alternatively  be  defined  using  the  wavefunction  pairs  (x^,  7d3)) 
and  (cr^\a^)  (Cormier,  1980).  This  fundamental  matrix  has  exactly  the  same  form  as  AI.4, 
but  with  (3)  replacing  the  (2)  superscripted  wavefunctions  and  the 'accent  replacing  the'  in  the 
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vertical  slownesses.  In  all  calculations,  the  (3)  superscripted  wavefunctions  are  substituted  for  the 
(2)  superscripted  (down-going)  wavefunctions  in  the  p  domains  in  which  exponentially  decaying 
and  growing  vertical  wavefunctions  exist.  With  a  few  simple  modifications  described  by  Doornbos 
(1981),  the  fundamental  matrix  defined  in  AI.4  can  be  applied  to  layers  having  a  negative  as  well 
as  a  positive  gradient  with  depth. 


Fundamental  Matrix  for  SH  Propagation 
The  SH  fundamental  matrix  and  its  inverse  are 


s; 

H 

t. 

fi~  1/2  ff(2) 

i  fi1/*  T)  (rf1) 

—  t  fj}/2  fff2) 

F(r)-1  =  V 
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Model  Parameterization 


(AI.5) 


Since  the  Langer  approximation  allows  layers  to  be  vertically  inhomogeneous,  the  effects  of  Earth 
curvature  are  built  into  the  model  parameterization.  All  formulae  are  evaluated  using  velocities 
and  densities  given  as  functions  of  radius,  r,  from  the  Earth’s  center.  In  each  inhomogeneous  layer, 
the  velocities  are  specified  by  analytic  functions,  which  have  only  one  turning  point  solution  in  the 
p  domain  of  interest.  Layer  boundaries  are  introduced  and  boundary  conditions  are  evaluated  at 
discontinuities  in  velocity  derivatives  as  well  as  first  order  discontinuities. 
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To  provide  analytic  forms  for  the  delay  time  functions  ra  and  Tp,  each  inhomogeneous  layer  is 
parameterized  by  making  the  flattened  velocity  be  a  linear  function  in  the  flattened  depth  coordi¬ 
nate,  z.  The  usual  (Muller,  1971)  mapping  between  the  flattened  velocity  function  Vf(z)  and  the 
true  velocity  function  v(r)  is  assumed  : 

v(r )  =  rvj(z)/Re 

where 

z~  =  Re  log(rc/ Re) 
where  Re  is  the  radius  of  the  Earth. 

The  flattened  velocity  function  vj  is  assumed  to  be  a  linear  function  in  flattened  depth,  com¬ 
puted  from  the  values  of  vj  at  flattened  depths  z~  and  z+_x  corresponding  to  radii  r~  and  r*_j, 
bounding  the  top  and  bottom,  respectively  of  vertically  inhomogeneous  layer  n.  The  analytic  form 
of  the  delay  time  function  r(r)  becomes 


r(r)  = 


zt  -  zl 


u/(*n)  “  ”/(*») 


IRIIV1  -  vj  -  Re/p 


In 


Relp  +  \JR-Hp1  -  vj' 


vf 


(AI.6) 


This  parameterization  is  adequate  in  representing  thick  regions  of  the  crust  and  uppermost 
mantle,  in  which  velocity  gradients  are  nearly  constant  or  slowly  varying.  Usually  less  than  ten 
inhomogeneous  layers  are  all  that  are  needed  to  describe  models  having  several  first  order  discon¬ 
tinuities  and/or  discontinuities  in  gradient.  Alternative  model  parameterizations,  which  give  an 
analytic  form  of  r,  are  discussed  by  Cormier  (1980),  Cerveny  and  Jansky  (1983),  and  Cormier  and 
Richards  ( 1989). 
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APPENDIX  II  -  Mode  Amplitudes  and  Eigenfunctions 


Rayleigh  Modes 

The  summation  of  locked  Rayleigh  modes  requires  the  calculation  of  an  antisymmetric  Y  matrix 
having  live  independent  elements. 

The  Y  Matrix 

At  the  radius  rc  at  the  top  of  the  capping  layer,  starting  values  of  the  Y  matrix  are  taken  to  be 

Yia  =  -  Al  -  B]  \ac 

Yi3  =  -Ac  p/re  -  Be  Xac  A/?c 

Yu  =  -i  pc  Xfc  (AII.l) 

Y23  =  i  pc  Xac 
Ym  =  -  Xp,  Xac  -  p7/rl 

where  i  =  V-T,  and 

=  *  \/p2/r?  “  l//?e 

=  *  'Jf/rl  -  l/o? 

Ac  =  2p2/r2  pc  -  pe 
Bc  =2  p2/r2  pc 
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and  pc,  Pei  ac,  (3C  are  the  shear  modulus,  density,  P  velocity,  and  S  velocity,  respectively,  of  the 
high  velocity  capping  layer. 

At  any  radius  r,  Y  can  be  computed  from  the  product 

Y(r)  =  KT(r,r+)Y(r;)K(r,r+)  (AII.2) 

where  K  is  a  P-SV  propagator  matrix  equal  to  a  product  of  intralayer  propagator  matrices  for  each 
layer,  m,  m  +  1,  etc. 

K  =  Km(r,r+)  Km+l(r~,  r++1) - Kn(r^_1,r+)  (AII.3) 

Layers  are  separated  by  boundaries  at  which  velocities  and/or  densities  have  either  first  or 
second  order  discontinuities.  Within  each  layer,  the  velocity  functions  are  continuous,  analytic 
functions.  Each  interlayer  propagator  matrix,  Km  is  constructed  from  the  zeroth  order  term  in 
frequency  of  the  uniform  asymptotic  approximation  to  the  fundamental  matrix  F  of  the  inhomo¬ 
geneous  layer.  Since  the  uniform  asymptotic  approximation  of  Langer  is  assumed,  the  velocity 
functions  within  each  layer  must  have  no  more  than  one  turning  point  for  each  ray  parameter,  p. 
With  this  restriction,  computations  can  still  be  conducted  in  a  complicated  model  having  one  or 
more  low  velocity  zones,  as  long  as  this  model  is  built  from  "layers”  in  which  the  analytic  functions 
for  P  and  S  velocity  have  only  a  single  turning  point  for  each  p. 

The  intralayer  propagator  is  defined  by 

Km(r-,r++1)  =  F(r- )  F"1(r++1)  (AII.4) 

Substituting  in  equation  AII.2  the  forms  for  the  fundamental  matrix  and  its  inverse  from 
equation  AI.4,  and  simplifying  the  resulting  expression  gives  recursion  relations  as  follows  for  the 
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upward  propagation  of  Y  matrix  elements: 


yi2(r")  =  *<*i(rn)  kWn  kGn  -  2  A(r+_,)  2?(>\t_i)  0Wn  o Gn 

*xl 

^i3(r;)  =  -  £  *<fc(r«)  [A(r+_1)S(r+_1)  -  p/r+.j  5(r„_i+)]  *TYn  *G„ 
*=i 

Yi^^n1)  =  -  td3(r~)  kWn  kGn 

k=  1 

*=1 

y34(r;)  =  -£  fcds(r“1)  *VFn  *Gn  +  2  p/r^i  o^n  oG„ 

k=l 

where  the  quantities  td/(z),  tWn,  kGn  are  defined  as  follows: 

*di(r)  =  -A2(r)  -  B2(r)  jfcA0(r) 
jkd2(r)  =  A(r)  p/r  -f-  5(z)  *Aa(r)  *A,j(r) 
kd3(r)  =  ipir)k\0(r) 
kd4{r)  =  -ip(r)  * AQ (r) 

*d5(r)  =*  A^(r)  *Aa(r)  +  (p/r)2 


Aiyn  =  -(kdi(r+_1)F34(r"_1)  +  2id2(rJ_1)Fi3(  ^_j)  +  fcd3(^_i)^i4(^-i 

+  fcd4(rit-l)^23(rn-l)  +  fc^5(rn-l)^12(rn-l) 


for  A:  0  and 


0VKn=  p/r+F12(r 


+  Kr^-i)  +  P/r»-i  t )]lri3(rn_ i ) 
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(All. 8) 


30 


V'2^r"Mrn-l) 

iGn  =  *(2)(0  a^(r-)  <r«(rti)  , 

yfipteWU) 

2 Gn  =  _,r(2)(rn  )  <7'(1)(rn  )  (AII.9) 

y/Mrn  )p(rU) 

3 Gn  =  -*(1)(r-)  a(2)(r")  ^(r+.r)  a(1)(r+_1) 

VWnM^-l) 

*Gn  =  ^(l)(rn)  ^(1)(r^)  ^(^-l)  ff(2>(^_l)  -7==^== 

VWW^t.  i) 

*Aa  and  denote  the  following  at  boundaries  r~  and  r+  : 
iA*(rn)  =  f(r“) 

iMtf-i)  =  tf'ti) 

lA/3(rn)  =  *?(0 
lM^-i)  =  77(r+_t) 

aM^n)  =  fc1) 

2A*(r+_1)  =  e(r+_x) 

2^p(r~)  =  ~v(rnl) 

2A0(r+_1)  =  ->?«_!) 

3Mrn)  =  -f(r«) 

3Aa(r+_i)  =  -f(rti) 

3^0(r~)  =  T}(r~) 

3^p(rn-l)  =  *Krn-l) 

<Mrn)  =  “?(»*“) 
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Mrt-i)  - 

4^(^n)  =  -»Krn) 
^('n-l)  =  “Ww-l) 


Layer  Reduction 

The  first  term  (Jb  =  1)  in  the  summation  in  equation  AII.5  is  of  the  same  form  as  the  starting 
values  Y  matrix  in  the  capping  layer  in  regions  of  slowness  in  which  the  vertical  wavefunctions 
behave  exponentially.  When  this  first  term  is  exponentially  larger  by  several  orders  of  magnitude 
than  the  ( k  —  2, 3, 4, 5)  terms,  then  the  Y  matrix  calculation  may  be  started  at  a  higher  layer, 
taking  this  higher  layer  as  the  capping  layer.  This  procedure  of  layer  reduction  is  analogous  to  that 
described  in  homogeneously  layered  models  (Panza  and  Sudhadolc,  1987). 

Eigenfunctions 

Although  propagation  of  the  Y  matrix  elements  has  been  shown  to  be  numerically  stable  at  ar¬ 
bitrarily  high  frequency  (Abo-Zena  1979;  Harvey,  1981),  numerical  problems  in  the  calculation  of 
the  Rayleigh  eigenfunctions  reoccur  if  E  is  calculated  by  multiplying  propagator  matrices.  One 
approach  to  this  problem  is  to  divide  a  layer  into  thin,  pseudo  layers,  and  rescale  the  propagator 
matrix  after  propagation  through  each  thin  layer.  Better  techniques,  however,  can  be  formulated, 
which  do  not  require  the  introduction  of  additional  pseudo  layers. 

One  technique,  described  by  Harvey  (1985),  expresses  the  eigenfunctions  in  terms  of  Y  matrix 
elements  by  propagating  the  wavefield  upward  from  the  cap  layer  and  downward  from  the  free 
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surface.  Thus,  since  the  calculation  of  Y  elements  is  numerically  stable,  so  is  the  calculation  of  the 
E  eigenfunctions.  In  this  technique,  eigenvalues  can  be  normalized  at  the  source  depth,  offering 
numerical  advantages  in  the  calculation  of  channel  waves  having  vanishingly  small  energy  outside 
of  a  waveguide. 

The  technique  used  here  also  does  not  require  pseudo  layering,  but  retains  the  standard  nor¬ 
malization  of  the  Ei  function  to  1  at  the  free  surface.  The  first  step  in  this  technique  is  to  recognize 
that  the  stress  eigenfunctions  £3  and  £4  can  be  calculated  from  the  displacement  eigenfunctions 
Ei  and  £2  by 


E3  =  -  Y14/Y34  Ei  -  Y24/Y34  £2 

(AII.10) 

E4  —  —  Y13/Y34  Ei  +  Y23/Y34  £2 

Using  these  relations,  the  four  equations  that  propagate  the  E  vector, 

E(r)  =  K(r,rn)  E(r„)  (AII.ll) 

can  be  rewritten  as  two  equations  that  propagate  £1  and  £2, 


and  the  two  equations  given  in 
functions. 


E\  (r) 

£2(r) 


-  L(z,zn ) 


Ei  (rn) 
Ei(  2) 


(AII.12) 


AII.10  between  the  displacement  eigenfunctions  and  stress  eigen- 
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A  new  2x2  propagator  matrix  L  is  defined  having  components 


L\ i  =  An  —  An  Yu/I'm  +  A14Y13/Y34 

In  =  Au  -  A13Y24/Y34  +  A14Y23/Y34 

(AII.13) 

£21  =  A21  —  A23  Y14/Y34  +  A24  Y13/Y34 

L22  =  A22  +  A23  Y24/Y34  +  A24  Y23/Y34 


To  ensure  numerical  precision  in  a  machine  calculation,  the  individual  propagator  elements  as 
well  as  the  recursion  formulae  in  All  5  for  the  Y  matrix  elements  must  be  substituted  into  the 
definitions  of  the  the  £<,  elements  in  AII.13,  a  fraction  formed  with  the  common  denominator  of 
Y34,  and  the  numerator  of  the  fraction  simplified.  When  this  simplification  is  done,  it  is  seen  that 
all  numerator  terms  that  potentially  are  of  the  largest  exponential  order  cancel.  Although  many 
cancellations  occur,  the  resulting  expressions  for  the  £tJ  elements  are  still  quite  lengthy  and  are 
not  given  here. 


Love  Modes 
D\  and  £>3 

In  this  case,  calculation  of  the  dispersion  function  D\  eigenfunction  vector  E  can  proceed  by  simple 
multiplication  of  propagator  matrices  without  loss  of  numerical  precision.  The  vector  {D\,  £>2 ) 
in  the  notation  of  Harvey  (1985)  is  equal  to  the  vector  Esh  in  the  notation  of  Cormier  (19S0). 
In  the  capping  layer,  ( D\ ,  £>2  is  simply  equal  to  the  first  row  of  the  inverse  fundamental  matrix 
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for  SH  waves.  Any  constant  may  be  chosen  to  multiply  the  starting  value  of  (D\,  D 2),  since  this 
constant  will  cancel  in  the  definition  of  eigenfunctions  and  in  the  ratio  ^ appearing  and 
the  expression  for  the  total  response.  Starting  values  of  D\  and  D2  at  the  top  of  the  cap  layer  are 
thus  taken  as 


D\  =  —  i  pc  Pc  ^0C 


(AII.14) 


D2  =  -1/& 

D\  and  D2  are  propagated  upward  by  multiplication  of  SH  propagator  matrices.  Since  {D\,  D2) 
are  related  to  the  inverse  fundamental  matrix,  one  must  right  multiply  the  starting  values  by  the 
SH  propagator  matrix. 


Di(r) 

A(re) 

D2(r) 

D2{rc) 

K  (r,r+) 


Eigenfunctions 

Love  wave  eigenfunctions  are  defined  by 


(AII.15) 


Ei(r)  =  D2{t)/D2{R e) 


(AII.16) 


E2{t)  =  D\{r)/ D2(Re) 


In  the  residue  calculation,  scale  factors  can  be  applied  in  each  layer  and  discarded  during 
upward  propagation.  This  is  because  all  scale  factors  cancel  when  ratio  >s  formed.  In 

the  eigenfunction  calculation,  the  total  scale  factor  of  each  £>,  must  be  saved  in  order  to  properly 
describe  regions  of  exponential  decay  of  the  eigenfunction.  In  the  cases  where  E\  and  Ei  are 
exponentially  small,  the  depth  of  the  capping  layer  can  be  raised  and  calculations  started  at  a 
shallower  depth. 


36 


p  (gm-cm'  *) 

Figure  2:  Test  models  having  three  different  intensities  of  gradients  in  an  inhomogeneous  layer 
overlaying  a  homogeneous  halfspace.  Model  1  is  the  test  model  of  Spudich  and  Ascher  (1983) 
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Figure  3:  (a)  Love  and  (b)  Rayleigh  mode  dispersion  curves  calculated  in  Model  1  using  a  thin 
layered  representation  of  the  gradient  layer  (solid)  and  the  Langer  approximation  in  a  continuous 
representation  of  the  gradient  layer  (dashed). 
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Figure  4:  Comparison  of  synthetic  seismograms  calculated  in  Model  1  using  a  thin  layered  represen¬ 
tation  of  the  gradient  layer  (solid)  and  the  Langer  approximation  in  a  continuous  representation  of 
the  gradient  layer  (dashed).  The  source  is  a  point  double  couple  at  4.92  km.  depth,  corresponding 
to  a  vertically  dipping  strike  slip  fault,  striking  to  the  north,  observed  at  receivers  at  45°  azimuth.  A 
step  function  time  dependence  of  the  scalar  moment  is  assumed.  Shown  are  the  three  components 
of  particle  velocity.  The  effects  of  geometric  spreading  of  body  waves  have  been  approximately 
removed  by  multiplying  each  seismogram  by  range. 
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Figure  5:  Comparison  of  synthetic  seismograms  calculated  in  Model  3  using  a  thin  layered  repre¬ 
sentation  of  the  gradient  layer  and  the  Langer  approximation  in  a  continuous  representation  of  the 
gradient  layer.  The  result  of  the  discrete  method  is  shown  at  each  range.  The  lower  amplitude 
trace  labeled  DIF  is  the  difference  between  the  seismograms  calculated  by  the  two  different  param- 
eterizations,  (D)  discrete  thin  layered  and  (CL)  continuous  with  the  Langer  approximation,  i.e., 
DIF{t)  =  Si j(t)  -  ScL(t)-  An  approximate  correction  for  geometric  spreading  of  body  waves 
has  been  made. 


41 


o  o  o  o  o  o’ 

O  CM  W  (0  CO  O 


(ui*)  qjdaa 


42 
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Figure  7:  A  comparison  of  the  synthetic  for  the  vertical  component  of  the  fundamental  mode 
Rayleigh  wave  using  perturbation  theory  and  an  exact,  complex  mode  calculation.  At  each  range, 
the  results  of  the  exact  calculation  are  followed  by  the  differential  seismogram  obtained  by  subtract¬ 
ing  the  seismogram  calculated  by  perturbation  theory  from  the  seismogram  calculated  by  complex 
modes  and  eigenfunctions.  Each  trace  is  normalized  by  its  peak  amplitude,  indicated  by  the  number 
to  the  left  of  each  trace. 
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Model  MQ 


Figure  8:  A  simple  crust  and  upper  mantle  model  MH  composed  of  two  homogeneous  crustal  layers 
and  overlaying  a  homogeneous  mantle;  and  model  MG  having  a  single  crustal  gradient  layer  and 
mantle  gradient  layer. 
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Figure  9:  A  comparison  of  synthetics  in  model  MH  (above)  and  MG  (below)  computed  by  the 
Langer-locked  mode  method,  summing  the  first  10  Rayleigh  modes.  Shown  is  the  vertical  dis¬ 
placement  for  a  double  couple  point  source  at  30  km  depth.  The  orientation  of  the  double  couple 
corresponds  to  a  vertically  dipping  strike  slip  fault,  striking  to  the  north,  observed  at  an  azimuth 
°f  45° •  A  step  function  time  dependence  of  the  scalar  moment  is  assumed,  and  the  result  has  been 
convolved  with  a  short  period  WWSSN  instrument  response. 


Riducid  Tim*  (sic),  t  -  r/  7.8< 


Figure  10:  A  comparison  of  synthetics  in  model  MH  (above)  and  MG  (below)  computed  by  the 
Langer-locked  mode  methods,  summing  all  of  the  Rayleigh  modes  in  a  frequency  band  up  to  2  Hz. 
Source,  receivers,  and  instrument  are  described  in  figure  9. 
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University  Park 

Los  Angeles,  CA  90089-0741 

Prof.  Larry  J.  Ruff 

Department  of  Geological  Sciences 

1006  C.C.  Little  Building 

University  of  Michigan 

Ann  Arbor,  MI  48109-1063 

Dr.  R.B.  Tittmann 

Rockwell  International  Science  Center 
1049  Camino  Dos  Rios 

P.O.  Box  1085 

Thousand  Oaks,  CA  9 1 360 

Dr.  Richard  Sailor 

TASC  Inc. 

55  Walkers  Brook  Drive 

Reading,  MA  01867 

Dr.  Gregory  van  der  Vink 

IRIS,  Inc. 

1616  North  Fort  Myer  Drive 

Suite  1440 

Arlington,  VA  22209 

John  Sherwin 

Teledyne  Geotech 

3401  Shiloh  Road 

Garland,  TX  75041 

Professor  Daniel  Walker 

University  of  Hawaii 

Institute  of  Geophysics 

Honolulu,  HI  96822 
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William  R.  Walter 
Seismological  Laboratory 
University  of  Nevada 
Reno,  NV  89557 


Dr.  Gregory  Wojcik 
Weidlinger  Associates 
4410  El  Camino  Real 
Suite  110 

Los  Altos,  CA  94022 

Prof.  John  H.  Woodhouse 
Hoffman  Laboratory 
Harvard  University 
20  Oxford  St. 

Cambridge,  MA  02138 

Prof.  Francis  T.  Wu 
Department  of  Geological  Sciences 
State  University  of  New  York 
at  Binghamton 
Vestal,  NY  13901 

Dr.  Gregory  B.  Young 
ENSCO,  Inc. 

5400  Port  Royal  Road 
Springfield,  VA  22151-2388 
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GOVERNMENT 


Dr.  Ralph  Alewine  III 

DARPA/NMRO 

1400  Wilson  Boulevard 

Arlington,  VA  22209-2308 

Paul  Johnson 

ESS-4,  Mail  Stop  J979 

Los  Alamos  National  Laboratory 
Los  Alamos,  NM  87545 

Mr.  James  C.  Battis 

GL/LWH 

Hanscom  AFB,  MA  01731-5000 

Janet  Johnston 

GL/LWH 

Hanscom  AFB,  MA  01731-5000 

Dr.  Robert  Blandford 

DARPA/NMRO 

1400  Wilson  Boulevard 

Arlington,  VA  22209-2308 

Dr.  Katharine  Kadinsky-Cade 
GL/LWH 

Hanscom  AFB,  MA  01731-5000 

Eric  Chael 

Division  9241 

Sandia  Laboratory 

Albuquerque,  NM  87185 

Ms.  Ann  Kerr 

IGPP,  A-025 

Scripps  Institute  of  Oceanography 
University  of  California,  San  Diego 
La  Jolla,  CA  92093 

Dr.  John  J.  Cipar 

GL/LWH 

Hanscom  AFB,  MA  01731-5000 

Dr.  Max  Koontz 

US  Dept  of  Energy/DP  5 

Forrestal  Building 

1000  Independence  Avenue 
Washington,  DC  20585 

Mr.  Jeff  Duncan 

Office  of  Congressman  Markey 

2133  Rayburn  House  Bldg. 

Washington,  DC  20515 

Dr.  W.H.K.  Lee 

Office  of  Earthquakes,  Volcanoes, 
&  Engineering 

345  Middiefield  Road 

Menlo  Park,  CA  94025 

Dr.  Jack  Evemden 

USGS  -  Earthquake  Studies 

345  Middiefield  Road 

Menlo  Park,  CA  94025 

Dr.  William  Leith 

U.S.  Geological  Survey 

Mail  Stop  928 

Reston,  VA  22092 

Art  Frankel 
USGS 

922  National  Center 
Reston,  VA  22092 


Dr.  T.  Hanks 
USGS 

Nat'l  Earthquake  Research  Center 
345  Middiefield  Road 
Menlo  Park,  CA  94025 

Dr.  James  Hannon 

Lawrence  Livermore  Nat’l  Laboratory 
P.O.  Box  808 
Livermore,  CA  94550 


Dr.  Richard  Lewis 

Director,  Earthquake  Engineering  &  Geophysics 
U.S.  Army  Corps  of  Engineers 
Box  631 

Vicksburg,  MS  39180 

James  F.  Lewkowicz 
GL/LWH 

Hanscom  AFB,  MA  01731-5000 


Mr.  Alfred  Lieberman 
ACDA/VI-OA'State  Department  Bldg 
Room  5726 
320 -21st  Street,  NW 
Washington,  DC  20451 
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Stephen  Mangino 
GL/LWH 

Hanscom  AFB,  MA  01731-5000 


Dr.  Frank  F.  Pilotte 

HQ  AFTAQTT 

Patrick  AFB,  FL  32925-6001 


Dr.  Robert  Masse 
Box  25046,  Mail  Stop  967 
.  Denver  Federal  Center 
Denver,  CO  80225 


Art  McGarr 

U.S.  Geological  Survey,  MS-977 
345  Middlefield  Road 
Menlo  Park,  CA  94025 


Richard  Morrow 
ACDA/VI,  Room  5741 
320  21st  Street  N.W 
Washington,  DC  20451 


Dr.  Keith  K.  Nakanishi 
Lawrence  Livermore  National  Laboratory 
P.O.  Box  808,  L-205 
Livermore,  CA  94550 


Dr.  Carl  Newton 

Los  Alamos  National  Laboratory 

P.O.  Box  1663 

Mail  Stop  C335,  Group  ESS-3 
Los  Alamos,  NM  87545 

Dr.  Kenneth  H.  Olsen 

Los  Alamos  Scientific  Laboratory 

P.O.  Box  1663 

Mail  Stop  C335,  Group  ESS-3 
Los  Alamos,  NM  87545 

Howard  J.  Patton 

Lawrence  Livermore  National  Laboratory 
P.O.  Box  808,  L-205 
Livermore,  CA  94550 


Mr.  Chris  Paine 
Office  of  Senator  Kennedy 
.  SR  315 

United  States  Senate 
Washington,  DC  20510 

Colonel  Jerry  J.  Perrizo 
AFOSR/NP,  Building  410 
Bolling  AFB 

Washington,  DC  20332-6448 


Katie  Poley 
CLA-OSWR/NED 
Washington,  DC  20505 


Mr.  Jack  Rachlin 
U.S.  Geological  Survey 
Geology,  Rm  3  Cl 36 
Mail  Stop  928  National  Center 
Reston.VA  22092 

Dr.  Robert  Reinke 
WL/NTESG 

Kirtland  AFB,  NM  87117-6008 


Dr.  Byron  Ristvet 

HQ  DNA,  Nevada  Operations  Office 
Attn:  NVCG 
P.O.  Box  98539 
Us  Vegas,  NV  89193 

Dr.  George  Rothe 
HQAFTAC/TGR 
Patrick  AFB,  FL  32925-6001 


Dr.  Alan  S.  Ryall,  Jr. 
DARPA/NMRO 
1400  Wilson  Boulevard 
Arlington,  VA  22209-2308 


Dr.  Michael  Shore 
Defense  Nuclear  Agency/SPSS 
6801  Telegraph  Road 
Alexandria,  VA  22310 


Dr.  Albert  Smith 

Los  Alamos  National  Uboratory 

L-205 

P.  O.  Box  808 
Livermore,  CA  94550 

Donald  L.  Springer 

Uwrence  Livermore  National  Uboratory 
L-205 

P.  O.  Box  808 
Livermore ,  CA  94550 
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Mr.  Charles  L.  Taylor 
GL/LWG 

Hanscom  AFB,  MA  01731-5(XX) 


DARPA/RMO/Security  Office 
1400  Wilson  Boulevard 
Arlington,  VA  22209 


Mr.  Steven  R.  Taylor 

Lawrence  Livermore  National  Laboratory 

L-205 

P.  O.  Box  808 
Livermore,  CA  94550 

Dr.  Eileen  Vergino 

Lawrence  Livermore  National  Laboratory 
L-205 

P.  O.  Box  808 
Livermore,  CA  94550 

Dr.  Thomas  Weaver 
Los  Alamos  National  Laboratory 
P.O.  Box  1663,  Mail  Stop  C335 
Los  Alamos,  NM  87545 


J.J.  Zucca 

Lawrence  Livermore  National  Laboratory 
P.  O.  Box  808 
Livermore,  CA  94550 


GL/SULL 
Research  Library 

Hanscom  AFB  ,  MA  01731-5000  (2  copies) 


Secretary  of  the  Air  Force 
(SAFRD) 

Washington,  DC  20330 


Office  of  the  Secretary  Defense 
DDR&E 

Washington,  DC  20330 


Geophysics  Laboratory 
Attn:  XO 

Hanscorr.  AFB,  MA  01731-5000 


Geophysics  Laboratory 
Attn:  LW 

Hanscom  AFB,  MA  01731-5000 


DARPA/PM 

1400  Wilson  Boulevard 

Arlington,  VA  22209 


Defense  Technical  Information  Center 
Cameron  Station 

Alexandria,  VA  22314  (5  copies) 


Defense  Intelligence  Agency 
Directorate  for  Scientific 
&  Technical  IntelligenceAttn:  DT1B 
Washington,  DC  20340-6158 


AFTAC/CA 

(STINFO) 

Patrick  AFB,  FL  32925-6001 


TACTEC 

Battelle  Memorial  Institute 
505  King  Avenue 

Columbus,  OH  43201  (Final  Report  Only) 


HQDNA 

Attn:  Technical  Library 
Washington,  DC  20305 


DARPA/RMO/RETREEVAL 
1400  Wilson  Boulevard 
Arlington,  VA  22209 
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rr>NT R ACTOR S  (Foreign) 


Dr.  Ramon  Cabre,  S.J. 
Observatorio  San  Calixto 
Casilla  5939 
La  Paz,  Bolivia 


Prof.  Hans-Peter  Harjes 
Institute  for  Geophysik 
Ruhr  University/Bochum 
P.O.  Box  102148 
4630  Bochum  1,  FRG 

ft 

Prof.  Eystein  Husebye 
NTNF/NORSAR 
P.O.  Box  51 

N-2007  Kjeller,  NORWAY 


Prof.  Brian  L.N.  Kennett 
Research  School  of  Earth  Sciences 
Institute  of  Advanced  Studies 
G.P.O.  Box  4 

Canberra  2601,  AUSTRALIA 

Dr.  Bernard  Massinon 

Societe  Radiomana 

27  rue  Claude  Bernard 

75005  Paris,  FRANCE  (2  Copies) 


Dr.  Pierre  Mecheler 
Societe  Radiomana 
27  rue  Claude  Bernard 
75005  Paris,  FRANCE 


Dr.  Svein  Mykkeltveit 
NTNF/inORSAR 
P.O.  Box  51 

N-2007  Kjeller,  NORWAY 
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FOREIGN  (Others) 


Dr.  Peter  Basham 

Earth  Physics  Branch 

Geological  Survey  of  Canada 

1  Observatory  Crescent 

Ottawa,  Ontario,  CANADA  K1A  0Y3 

Dr.  Eduard  Berg 
Institute  of  Geophysics 
University  of  Hawaii 
Honolulu,  HI  96822 


Dr.  Michel  Bouchon 
I.R.l.G.M.-B.P.  68 
38402  St.  Martin  D'Heres 
Cedex,  FRANCE 


Dr.  Hilmar  Bungum 
NTNF/NORSAR 
P.O.  Box  5 1 

N-2007  KjeUer,  NORWAY 


Dr.  Michel  Campillo 
Observatoire  de  Grenoble 
I.R.1.G.M.-B.P.  53 
38041  Grenoble.  FRANCE 


Dr.  Kin  Yip  Chun 
Geophysics  Division 
Physics  Department 
University  of  Toronto 
Ontario,  CANADA  M5S  1A7 

Dr.  Alan  Douglas 
Ministry  of  Defense 
Blaclcncst,  Brimpton 

Reading  RG7-4RS,  UNITED  KINGDOM 


Dr.  Fekadu  Kebede 
Seismological  Section 
Box  12019 

S-750  Uppsala,  SWEDEN 


Dr.  Tormod  Kvaema 
NTNF/NORSAR 
P.O.  Box  51 

N-2007  KjeUer,  NORWAY 


Dr.  Peter  Marshal 
Procurement  Executive 
Ministry  of  Defense 
Blacknest,  Brimpton 

Reading  FG7-4RS,  UNITED  KINGDOM 

Prof.  An  Ben-Menahem 
Department  of  Applied  Mathematics 
Weizman  Institute  of  Science 
Rehovot,  ISRAEL  951729 


Dr.  Roben  North 

Geophysics  Division 

Geological  Survey  of  Canada 

1  Observatory  Crescent 

Ottawa,  Ontario,  CANADA  K1A0Y3 

Dr.  Frode  Ringdal 
NTNF/NORSAR 
P.O.  Box  51 

N-2007  Kjeller,  NORWAY 


Dr.  Jorg  Schlittenhardt 

Federal  Institute  for  Geosciences  &  Nat’l  Res. 
Postfach  510153 

D-3000  Hannover  51,  FEDERAL  REPUBLIC  OF 
GERMANY 


Dr.  Roger  Hansen 
NTNF/NORSAR 
P.O.  Box  51 

N-2007  KjeUer.  NORWAY 


Dr.  Manfred  Henger 

Federal  Institute  for  Geosciences  &  Nat'l  Res. 
Postfach  510153 
D-3000  Hanover  51,  FRG 


Ms.  Eva  Johannisson 
Senior  Research  Officer 
National  Defense  Research  Inst. 

P.O.  Box  27322 

S-102  54  Stockholm,  SWEDEN  -12- 


